Monte Carlo Simulation of a Model of Water 
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We simulate TIP3P water using a constrained Monte Carlo algorithm to generate electrostatic interactions 
eliminating the need to sum over long ranged Coulomb interactions. We study discretization errors when in- 
terpolating charges using splines and Gaussians. We compare our implementation to molecular dynamics and 
Brownian dynamics codes. 
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The TIP3P model of water 1 1] is often used to study the 
accuracy of algorithms for atomistic simulation. The model 
has a single Lennard-Jones center representing an oxygen 
atom together with three charges (-0.834, +0.417, +0.417) 
arranged in a triangle. The oxygen-hydrogen bond length is 
0.9572A the angle between bonds is 104.52°. Accurate sim- 
ulation of this model is surprisingly challenging: The bare 
electrostatic interaction between oxygen atoms at a separation 
of 2. 75 A, is over 100 fcsT. Small errors in the representation 
of the electrostatic potential lead to significant errors in the 
total energy due to large cancellations; the binding energy per 
hydrogen bond is only 7 ksT. 

Many molecular dynamics codes for the simulation of large 
numbers of charges are based on Poisson solvers. The codes 
interpolate charges to a cubic grid and then calculate the elec- 
trostatic energy via fast Fourier transform |2| or multigrid 
fs, 4). The principle difficulty is controlling errors in the 
Coulomb interaction using high order interpolation. One re- 
quires a relative eiTor of at most ~ 10"'*. In this article we 
present a Monte Carlo algorithm for simulation at this level of 
accuracy. We avoid solving the Poisson equation by general- 
izing an algorithm which generates the Coulomb interaction 
between particles using Monte Carlo evolution of the elec- 
tric field. Previous codes using this local algorithm have been 
of low accuracy, sufficient for the study of lattice gasses ISJ 
or charges interacting through an implicit solvent |6]. They 
were still far from the accuracy needed for the simulation of 
TIP3P. This articles considers the modifications necessary to 
the algorithm in order to reliably simulate standard atomistic 
models. 

There were three important sources of error in the energy 
functions used in previous work with local electrostatics algo- 
rithms: 

• use of low order interpolation leading to distorted 
charge distributions 

• aliasing errors in the lattice Green functions leading to 
a self-energy with the periodicity of the lattice 

• low order discretization of the lattice Green function 
leading to anisotropy in the effective interactions. 

In many codes interpolation of charges from the continuum 
to the cubic grid is performed with splines. A one-dimension 



ri-spine is a set of n polynomials of order n — I. These poly- 
nomials give the quantity of charge which is deposited on n 
consecutive sites of the lattice as a function of the position 
of the particle, fi{x), 1 < i < n. Linear interpolation cor- 
responds to a 2-spline. In three dimensions one takes the 
product of splines in the x, y and z directions, interpolating 
a charge to lattice site, thus /i(r) = fi{x)fj{y)fk{z) for 
r = {x,y,z) and I = {i,j,k). Splines have several useful 
properties for interpolation: They conserve total charge ex- 
actly; they are smooth with n — 2 continuous derivatives. With 
Fourier solvers splines work well if one takes n > 4 |2]. 

An alternative to splines is interpolation with truncated 
Gaussians |2, 3]. Consider interpolating a unit charge to a 
one-dimensional grid with Gaussian interpolation: fi{x) — 
exp {—{x — iy /2a'^)/^/2TT(j. The total interpolated charge 
can be evaluated for a large with the Poisson resummation 
formula: qi^t = Y.i M^) = Y.pfi'^'^P) where / is the 
(continuous) Fourier transform of fi{x). We find qint ^ 
1 + 2 cos(27ra;)e"^'^ . Already for (7 = 1 eiTors in charge 
conservation are 0(10"*). In practice one truncates beyond 
Act where A ^ 4 — 5, leading to an additional error which 
varies as e 

In order to study the various errors generated with lattice 
Monte Carlo algorithms for the electrostatic energy consider 
the interaction between two particles placed at r and r'. 

C/(r,r') = ^/,(r)G(l-m)/^(r') (1) 

l,m 

= E/(§3/(q-2-P)/(q)G(q) 

^g*q-(r-r')+2,.»p-r 

G{1) is the lattice Green function of the interaction between 
two sites and G(q), its Fourier transform, has the periodicity 
of the Brillouin zone. 27rp is a vector of the reciprocal lattice. 
The integral is over all Fourier space. The particles also have 
a self energy C/(r, r)/2. 

Consider the contribution p = in eq. (|2j 

C/o(r,0) = J G(q)/2(q)e^<i--^ (3) 
If / is Gaussian and G{r) = l/iirr, G{q) = 1/q'^ so that 

U„ir.O).'-^,^ (4) 
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This is the central formula for so called "particle mesh Ewald" 
methods \2, 3]. One neglects contributions with p 7^ and 
calculates the Coulomb energy as a lattice energy, eq. Q, 
plus a short range correction. Deviations from eq. (|4} oc- 
cur if the structure factors are not Gaussian: For splines 
/ = Yla siiic"(qQ,/2) with a — {x, y, z) which has the cumu- 
lant expansion / ^ J| exp(— — nq^/2880). Splines 
converge for large n to Gaussians of width cr^ = n/12. How- 
ever interpolation with splines generates an extra contribution 
at q = in eq.Q due to the term in in the cumulant. This 
leads in real space to an error which decays as l/r^. The am- 
plitude of this error decays rather slowly with n. 

The aliasing error comes from the contributions p ^ 0. 
Consider, for instance the self energy U{r, r)/2 and the con- 
tribution to eq. (|2} from pi ~ (1,0,0). Since /(q) de- 
cays rapidly in Fourier space the product /(q — 27rpi)/(q) 
is maximum on the boundary of the first Brillouin zone near 
q — 7r(l,0,0). If we sum over all symmetry related lattice 
reciprocal vectors we find a periodic one body potential 



Vi - /^(7rpi)G'(7rpi)2_]cos27rra 



(5) 



Higher order corrections to Vi come from larger p. We 
compare spline and Gaussian interpolation: For a Gaussian 
/^(ttpi) — e^'^ , whereas for a n-spline we find = 
(2/7r)2" ~ e"°-^". Requiring p - lO"'^ impHes that - 1 
or n ^ 10. An implementation using low order splines with 
only n — showed strong aliasing artefacts (7|. The sinu- 
soidal form of eq. (|5} permits simple analytic subtraction, but 
we will not pursue this point here. 

We now turn to errors in the lattice Green function, G(r). 
Coulomb's law in electromagnetism results from the impo- 
sition of a linear constraint. Gauss' law: V E = p, on a 
quadratic energy functional: U ~ 1/2 J (Pr. Previous 
codes that discretized these equations led to the standard 7- 
point discretization of the Laplacian operator: 
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Expanding we find 



^(q) = 2^(i 



cosqa) 



(6) 
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The presence of terms which involve q^/q^ imply a correc- 
tion to G(r) which decays as only l/r^. We now construct 
a discretization which converges faster Consider an energy 
which is a quadratic function of P electric field variables Ei 
where the subscript includes both positional and directional 
information. 



1 ^ 



EKE 12 



(8) 



We continue with an operator notation for compactness. We 
submit this energy to L < P linear constraints, ci 



ci 



ipEp 



= 0, V/ 



(9) 



where Dip is a discretization of the divergence operator at I 
and ei is the charge. Stationary points are found by consider- 
ing the functional A = Ue — J2i 't'lCi- D \s2l linear operator, 
we define the adjoint D* . The variational equations of E are 



KE~ D* 







(10) 



Solving for E in eq. (I10> and substituting in the constraint 
equation eq.(|9} we find a generaUzed Poisson equation for the 
Lagrange multipliers 



DR-^D* 







(11) 



From this Poisson equation we find that the minimum of the 
constrained energy is 



= ecj)/2 = eGKe/2 



(12) 



with the Green function ~ DK^^D* . We make contact 
with electrostatics if we recognize that D = div implies D* = 
— grad, and G^^ — — V^. 

We will now generalize these results to non-zero tempera- 
tures and show that the effective interaction between particles 
is still described by the Green function Gk- The constraints 
are now imposed by delta-functions in a partition function 

. p L 
Z= XldEpe-P^'^llSici) (13) 

P=i 1=1 

We decompose the field E into generalized "longitudinal" and 
"transverse" components by writing E — K^^D*(p + Et and 
change integration variables from Eto Ef. The partition func- 
tion then factorizes 



Z - e-'^^^ 



/ f[dE,,pe-^^^^^^/^f[S{DE, 

p=i 1=1 



Zji X constant 



(14) 



Where is a partition function for particles interacting with 
the Green function Gk- 

Previous implementations 0| took the following forms for 
the operators D and K: DE was the flux out of the site I to the 
six nearest neighbor sites. K was diagonal, K — Sij, leading 
to eq. (|6|l. Here we keep the same form for the operator D 
but for K we include interactions between neighboring links 
on the lattice; for x-oriented bonds of the lattice the energy 
function is 



Ue 



_F2 



1 

12^ 



tE, 



'i.j,k\x {-^i-\-l,j.k\x 



^-^i—l,j.k\x ) 

(15) ^ 
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with similar expression for the Hnks in the y and z directions. 
In Fourier space we find that 

K{q) = idiag(5 + cosgj;, 5 + cosgj^, 5 + cosg^) 
D(q) = (l-e*«^l-e^'^y^ 
L»*(q) = (l-e-*«^l-e-'■'^^l-e-*'^.)T 

where "diag" denotes a matrix with the indicated diagonal 
elements, so that 

5 + cos Oa 

OL 

G(q) = ^ + AE^ + --- (16) 

This form of G leads to reduced artefacts in the lattice Green 
function; errors now decay as 




r/o 

FIG. 1 : Scaled error in the pair potential using the energy eq. <16> 
for cr = 0.90, 1.0, 1.2, 1.4, 1.6. One particle placed at (0, 0, 0) the 
second displaced in the direction (1, 1, 1). Solid lines: Curves for 
CT = 1.2, 1.4, 1.6 superpose: errors in G dominate. Also included 
a single curve for cr = 1.2 starting at (0, 0, 0.5) which scales in an 
identical manner. Dashed lines: a = 0.90, 1.0. Oscillations, eq. Js}, 
from aliasing are also important and violate the scaling in a^. 

To calibrate the effective interaction generated by our con- 
strained algorithm we numerically inverted the Green function 
eq. (I16> . We take two interpolated unit charges and measured 
the potential between them. Figure Q, as a function of a and 
compared with the (exact) Ewald energy, Uw We find col- 
lapse of the error when we plot 47r(C/ — Uw)cr^ as a function 
of r/a for a > 1. We conclude that the error in the pair po- 
tential can be written in the scaling form 

dUG{r,r',a)^^V5{{r^r')/a) (17) 

for a > 1. V5 does depend on the direction of the relative 
displacement with respect to the lattice. The error increases 
strongly for r/a < 2, however smaller distances in our sim- 
ulations will be within the core of the Lennard- Jones poten- 
tial and will not be sampled. Due to the regularity of V5 one 



could also improve accuracy of the simulation by subtracting 
the error off of the real space potential after parameterizing it 
with splines. For a < 1 aliasing eiTors are increasingly impor- 
tant, adding an sinusoidal contribution to the eiTor as expected 
from eq. (|5j- cr = 1 generates eiTors in the potential which are 
0(10-'*). 

One simulates a system described by the energy eq. (jSJl with 
the constraint eq. (|9} with the Metropolis method by introduc- 
ing two independent Monte Carlo updates: Plaquette updates, 
which satisfy D5E = 0, consist of the coupled update of the 
4 links forming a plaquette fs*] of the cubic lattice. On each 
of these links the field is modified by the same amount A, so 
that the flux of E at each node remains constant. With the 
addition of nearest neighbor interactions, eq. il5\ . calculation 
of the energy change requires the values of the field from 12 
links. Motion of a particle is possible if a local update of the 
field is performed simultaneously such that D5E = Sp where 
6p is the localized charge fluctuation. 

We implemented a simulation of TIP3P water using Gaus- 
sian interpolation due to the superior convergence properties 
at higher accuracies. We work in units of the mesh size. Three 
dimensional Gaussians, with a = 1, are calculated as di- 
rect products of one dimensional Gaussians each truncated at 
A = 4.3. The three atoms of each molecule are interpolated 
together. The small error, 0(10^'*), in charge conservation is 
coiTected on the grid point nearest the oxygen atom. In this 
way we insure that charge is conserved in the algorithm to 
machine precision. We perform a trial move and re-perform 
both interpolation and charge coiTection steps. This gives us 
a localized charge fluctuation Sp. We generate a local field 
modification in a box enclosing the original and final sites: 
We use a SE such that J^i ^E^f /2 is minimum while respect- 
ing D5E = 5p. Outside of the region where 5p ~ we 
impose that 5E = 0. This leads to a small Poisson prob- 
lem (with zero flux boundary conditions) within the interpola- 
tion volume which can be solved using the FFTW library j^l- 
Lennard- Jones and crfc interactions were truncated at a dis- 
tance of 9A. We use Monte Carlo updates in both the position 
and orientation of the molecules, tuned to give an acceptance 
rate of about 40%. For each update of a molecule we perform 
100 plaquette updates. Simulations were performed at 30QK 
at constant volume. Due to the simplicity of the plaquette up- 
dates compared with the the calculation of the erfc interaction, 
they take a small part of the CPU time. 

We compared our Monte Carlo simulation of TIP3P with a 
molecular dynamics simulation 1 10] using using a Langevin 
thermostat, friction 1 ps^^, integration time step 1 fs. We 
used a cubic box of side 18.62 A and a grid of 20^ sites for 
the Monte Carlo. We measured the autocorrelation time of 
the potential energy, V using blocking |11], after removing 
the energy in the transverse electric field in the Monte Carlo 
runs. A set of recordings corresponding to T sweeps or time 
steps is averaged in blocks of 6 = 2™ recordings, to esti- 
mate the mean potential energy, {V), and a running estimate 
in the error in {V), a{h). For large blocking factors a{b) sat- 
urates to a constant a^', the integrated autocorrelation time is 
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given by r = alT/2{V'^ — (^)^)- In order to obtain reason- 
able statistics for the dynamics we simulated a small system 
of 216 particles for several thousand r. The running estimate 
of t(6) is plotted in Figure|2] We estimate that t = 1100 for 
molecular dynamics and t = 800 for Monte Carlo. We also 
performed Brownian dynamics simulations with the time step 
equal to 1/10 of the stability limit using an Euler integrator 
finding r = 3200. 
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FIG. 2: Blocking analysis of the energy with Monte Carlo (MC), 
molecular dynamics (MD) and Brownian dynamics (BD) to estimate 
the energy autocorrelation time. The estimated r for each method is 
indicated with a dot-dashed line. 

To conclude, we have implemented a Monte Carlo algo- 
rithm for the simulation of TIP3P. Each Monte Carlo time 
step implies two interpolations per charge plus calculation of 
a localized current. Each molecular dynamics step requires 
one charge interpolation and then three extrapolations for the 
force plus solution of the Poisson equation. The total com- 
plexity of the interpolation steps is very similar in molecu- 
lar dynamics (particularly multigrid) codes and in our Monte 
Carlo formulation. The integrated autocorrelation time with 
our algorithm is comparable to simulations performed using 
molecular dynamics when measured in sweeps. CPU time 
comparisons were less favorable to our code since the (Fourier 
based) Gromacs package contains contains extensive assem- 
bly routines for interpolation which we did not implement; the 
present Monte Carlo implementation is an order of magnitude 
slower. 

The authors in J3| reduced the range of the interpolation 



step by introducing a second on-lattice convolution after in- 
terpolating charges to the grid. Similar techniques are possi- 
ble here if we modify the kernel K so as to include a extra 
spreading step; charge interpolation then becomes cheaper. In 
our algorithm this simplification in charge motion is balanced 
by an increase in the complexity of plaquette updates. 

Our algorithm has the important advantage over other codes 
of being purely local, and thus easily implemented on parallel 
computers with limited interprocessor communication, such 
as is the case on low cost clusters. With a Monte Carlo algo- 
rithm there is also an enormous gain in flexibility in heteroge- 
nous environments: In simulations of a biomolecule or an in- 
terface most of the water molecules play the role of distant 
spectator, even if they provide the majority of charge centers. 
In Monte Carlo it is trivial to bias moves towards interesting 
degrees of freedom, or even introduce cluster or multi- 
step fl3"] updates; with molecular dynamics multi-scale and 
multi-step algorithms are difficult to implement and prone to 
instability. 
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